#### read in the 30 percent data
require(dplyr)
require(tidyr)
require(stargazer)
require(rsq)

fol_path= "D:\\Research\\global_refor_mac\\tabulated_data\\full_data\\"

# read in datasets
df_afr = readRDS(paste0(fol_path, "afr_30per_full.rds"))
df_asi = readRDS(paste0(fol_path, "asi_30per_full.rds"))
df_lat = readRDS(paste0(fol_path, "lat_30per_full.rds"))

df_plant = readRDS(paste0(fol_path, "plantat_long_agg.rds"))

# create new layers of forest cover
df_afr$per_forest_2000 = df_afr$per_all_forest + df_afr$per_for_loss
df_afr$per_forest_2000_2 = df_afr$per_forest_2000^2
df_afr$per_forest_2000_3 = df_afr$per_forest_2000^3
df_afr$per_forest_2000_4 = df_afr$per_forest_2000^4

# calculate percent 
df_asi$per_forest_2000 = df_asi$per_all_forest + df_asi$per_for_loss
df_asi$per_forest_2000_2 = df_asi$per_forest_2000^2
df_asi$per_forest_2000_3 = df_asi$per_forest_2000^3
df_asi$per_forest_2000_4 = df_asi$per_forest_2000^4

df_lat$per_forest_2000 = df_lat$per_all_forest + df_lat$per_for_loss
df_lat$per_forest_2000_2 = df_lat$per_forest_2000^2
df_lat$per_forest_2000_3 = df_lat$per_forest_2000^3
df_lat$per_forest_2000_4 = df_lat$per_forest_2000^4

# create left join of data
df_lat = left_join(df_lat, df_plant, by = "uni_ID")
df_afr = left_join(df_afr, df_plant, by = "uni_ID")
df_asi = left_join(df_asi, df_plant, by = "uni_ID")

# create percentage
df_lat$per_shc_all = df_lat$per_shc_1 + df_lat$per_shc_2 + df_lat$per_shc_3 + df_lat$per_shc_4 
df_afr$per_shc_all = df_afr$per_shc_1 + df_afr$per_shc_2 + df_afr$per_shc_3 + df_afr$per_shc_4 
df_asi$per_shc_all = df_asi$per_shc_1 + df_asi$per_shc_2 + df_asi$per_shc_3 + df_asi$per_shc_4 

df_eco = readRDS(paste0(fol_path, "eco_region.rds"))
csv_eco_region = read.csv(paste0("D:/Research/global_refor_mac/data/ecoregion_data/",
                                 "ecoregion_biome_translation.csv"), stringsAsFactors = FALSE)
df_eco$biome_num = csv_eco_region$BIOME_NUM[match(df_eco$ecoregions, csv_eco_region$ECO_ID)]

df_lat = left_join(df_lat, df_eco, by = "uni_ID")
df_afr = left_join(df_afr, df_eco, by = "uni_ID")
df_asi = left_join(df_asi, df_eco, by = "uni_ID")

df_lat$cont_var = "LatinAmerica"
df_afr$cont_var = "Africa"
df_asi$cont_var = "Asia"

# create full data
df_full = bind_rows(df_lat, df_afr, df_asi)

df_full$biome_num_f = factor(df_full$biome_num)
df_full$cont_var_f = factor(df_full$cont_var)
df_full$cont_num_f = factor(ifelse(df_full$cont_var == "LatinAmerica", 1,
                            ifelse(df_full$cont_var == "Africa", 2, 3)))
df_full$ecoregions_f = factor(df_full$ecoregions)
df_full$country_f = factor(df_full$country)

# save out the rds files 
# saveRDS(df_full, paste0("D:/Research/global_refor_mac/results/rds_reg_files/", "all_cont_reg.rds"))

# run continental regression 
reg_loss_cont = rxGlm(per_for_loss ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                      per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                      elev + iucnlow_2000 + iucnhigh_2000 + F(biome_num) + F(cont_num_f), 
                    data = df_full[(df_full$per_forest_2000 < 1), ], 
                    family = poisson(link = log), 
                    coefLabelStyle = "R", 
                    maxIterations = 1000)

pred_reg = rxPredict(reg_loss_cont, df_full, type = c("response"))

reg_loss_cont_glm <- glm(per_for_loss ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                    per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                    elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f + cont_num_f, 
                  family="poisson", 
                  data=df_full[df_full$per_forest_2000 < 1, ])


reg_gain_cont = rxGlm(per_for_gain ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                        per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                        elev + iucnlow_2000 + iucnhigh_2000  + F(biome_num) + F(cont_num_f), 
                      data = df_full[df_full$per_forest_2000 > 0, ], 
                      family = poisson(link = log), 
                      coefLabelStyle = "R", 
                      maxIterations = 1000)

reg_gain_cont_glm <- glm(per_for_gain ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                           per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                           elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f + cont_num_f, 
                         family="poisson", 
                         data=df_full[df_full$per_forest_2000 > 0, ])

# save out the rds files 
saveRDS(reg_loss_cont, paste0("D:/Research/global_refor_mac/results/rds_reg_files/", "reg_loss_cont_v1.rds"))
saveRDS(reg_gain_cont, paste0("D:/Research/global_refor_mac/results/rds_reg_files/", "reg_gain_cont_v1.rds"))


##########################################################################
### create stargazer results 
library(stargazer)
out_path = "D:/Research/global_refor_mac/results/regression_html/"

rsq_vec = c(rsq(reg_loss_cont_glm), rsq(reg_gain_cont_glm))
rsq_vec = round(rsq_vec, digits = 2)

stargazer(reg_loss_cont_glm, reg_gain_cont_glm,
          type = "html",
          title = "Global land use change model of deforestation and reforestation",
          out = paste0(out_path, "global_land_use_cont_biome.html"),
          p.auto = TRUE,
          align=FALSE, 
          digits = 5,
          font.size = "small",
          star.cutoffs=c(0.10,0.05,0.01),
          dep.var.labels=c("Percent deforestation", "Percent reforestation"),
          model.numbers = TRUE,
          multicolumn = FALSE,  
          perl = TRUE,
          keep = c(names(reg_loss_cont_glm$coefficients)[c(2:11, 23,24,1)]),
          covariate.labels= c("Agricultural revenue ( hec)",
                              "Percent forest cover",
                              "Percent forest cover 2",
                              "Percent forest cover 3",
                              "Percent forest cover 4",
                              "Distance to city (km)",
                              "Slope (degree)",
                              "Elevation (m)",
                              "Multiple protected area (percent of cell)",
                              "Strict protected area (percent of cell)",
                              "Africa FE",
                              "Asia FE",
                              "Intercept"),
          no.space = TRUE,
          keep.stat=c("n"),
          single.row = FALSE,
          add.lines = list(c("R-squared", rsq_vec),
                           c("Biome FE", "Yes", "Yes")))

######################################################################
### create hardwood regression 

hardwood_exp = read.csv("D:/Research/global_refor_mac/data/hardwood/hardwood_export.csv",
                        stringsAsFactors = FALSE)
colnames(hardwood_exp) = c("export_value", "country")

hard_df = left_join(df_full, hardwood_exp, by = c("country" = "country"))


# run continental regression 
reg_loss_cont = rxGlm(per_for_loss ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                        per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                        elev + iucnlow_2000 + iucnhigh_2000 + export_value + F(biome_num) + F(cont_num_f), 
                      data = hard_df[(hard_df$per_forest_2000 < 1), ], 
                      family = poisson(link = log), 
                      coefLabelStyle = "R", 
                      maxIterations = 1000)

reg_gain_cont = rxGlm(per_for_gain ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                        per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                        elev + iucnlow_2000 + iucnhigh_2000 + export_value  + F(biome_num) + F(cont_num_f), 
                      data = hard_df[hard_df$per_forest_2000 > 0, ], 
                      family = poisson(link = log), 
                      coefLabelStyle = "R", 
                      maxIterations = 1000)


# save out the rds files 
saveRDS(reg_loss_cont, paste0("D:/Research/global_refor_mac/results/rds_reg_files/", "reg_loss_hardwood.rds"))
saveRDS(reg_gain_cont, paste0("D:/Research/global_refor_mac/results/rds_reg_files/", "reg_gain_hardwood.rds"))

######################################################################
## run the regressions

# run continental regression 
reg_loss_cont = rxGlm(per_for_loss ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                        per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                        elev + iucnlow_2000 + iucnhigh_2000 + F(country_f) + F(ecoregions_f), 
                      data = df_full[(df_full$per_forest_2000 < 1), ], 
                      family = poisson(link = log), 
                      coefLabelStyle = "R", 
                      maxIterations = 1000)

reg_gain_cont = rxGlm(per_for_gain ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                        per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                        elev + iucnlow_2000 + iucnhigh_2000 + F(country_f) + F(ecoregions_f), 
                      data = df_full[(df_full$per_forest_2000 > 0), ], 
                      family = poisson(link = log), 
                      coefLabelStyle = "R", 
                      maxIterations = 1000)

# save out the csv file
fact_df = data.frame("country" = levels(df_full$country_f),
                     "country_fixed" = paste0("F_country_f", 1:length(levels(df_full$country_f))),
                     stringsAsFactors = FALSE)

eco_df = data.frame("ecoregion" = levels(df_full$ecoregions_f),
                    "ecoregion_fixed" = paste0("F_ecoregions_f", 1:length(levels(df_full$ecoregions_f))),
                    stringsAsFactors = FALSE)



# save out the rds files 
saveRDS(reg_loss_cont, paste0("D:/Research/global_refor_mac/results/rds_reg_files/", "reg_loss_eco_country.rds"))
saveRDS(reg_gain_cont, paste0("D:/Research/global_refor_mac/results/rds_reg_files/", "reg_gain_eco_country.rds"))

saveRDS(eco_df, paste0("D:/Research/global_refor_mac/results/rds_reg_files/", "eco_fixed_eff.rds"))
saveRDS(fact_df, paste0("D:/Research/global_refor_mac/results/rds_reg_files/", "country_fixed_eff.rds"))


reg_loss_cont_glm <- glm(per_for_loss ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                           per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                           elev + iucnlow_2000 + iucnhigh_2000, 
                         family="poisson", 
                         data=df_full[df_full$per_forest_2000 < 1, ])

reg_gain_cont_glm <- glm(per_for_gain ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                           per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                           elev + iucnlow_2000 + iucnhigh_2000, 
                         family="poisson", 
                         data=df_full[df_full$per_forest_2000 > 0, ])

reg_gain_cont_glm$coefficients[1:11] = reg_gain_cont$coefficients[1:11]
std_gain = diag(vcov(reg_gain_cont_glm))
std_gain = reg_gain_cont$coef.std.error[1:11]
names(std_gain) = names(reg_gain_cont_glm$coefficients)

reg_loss_cont_glm$coefficients[1:11] = reg_loss_cont$coefficients[1:11]
std_loss = diag(vcov(reg_loss_cont_glm))
std_loss = reg_loss_cont$coef.std.error[1:11]
names(std_loss) = names(reg_loss_cont_glm$coefficients)



out_path = "D:/Research/global_refor_mac/results/regression_html/"

stargazer(reg_loss_cont_glm, reg_gain_cont_glm,
          type = "html",
          se = list(std_loss, std_gain),
          title = "Global land use change model of deforestation and reforestation",
          out = paste0(out_path, "global_land_use_cont_biome_count_ecoregion.html"),
          p.auto = TRUE,
          align=FALSE, 
          digits = 5,
          font.size = "small",
          star.cutoffs=c(0.10,0.05,0.01),
          dep.var.labels=c("Percent deforestation", "Percent reforestation"),
          model.numbers = TRUE,
          multicolumn = FALSE,  
          perl = TRUE,
          keep = c(names(reg_loss_cont_glm$coefficients)[c(2:11, 1)]),
          covariate.labels= c("Agricultural revenue ( hec)",
                              "Percent forest cover",
                              "Percent forest cover 2",
                              "Percent forest cover 3",
                              "Percent forest cover 4",
                              "Distance to city (km)",
                              "Slope (degree)",
                              "Elevation (m)",
                              "Multiple protected area (percent of cell)",
                              "Strict protected area (percent of cell)",
                              "Intercept"),
          no.space = TRUE,
          keep.stat=c("n", "rsq"),
          single.row = FALSE,
          add.lines = list(c("Country FE", "Yes", "Yes"),
                           c("Ecoregion FE", "Yes", "Yes")))



#################################################################
#### create regression that excludes deserts and mangroves

# read in biome category dataset
biome_cat = read.csv("D:/Research/global_refor_mac/data/carbon_file/biome_cat.csv",
                     stringsAsFactors = FALSE)
npv_carbon_refor = read.csv("D:/Research/global_refor_mac/data/carbon_file/npv_carbon_refor_trans.csv",
                            stringsAsFactors = FALSE)

# read in the carbon stock per decade
carbon_stock_dec = read.csv("D:/Research/global_refor_mac/data/carbon_file/carbon_stock_per_decade.csv",
                            stringsAsFactors = FALSE)
colnames(carbon_stock_dec)[4:13] = 1:10

# calculate biome reforestation category 
biome_three_var = biome_cat$biome_cat[match(df_full$biome_num,
                                            biome_cat$biome_num)]
# create continent variable
cont_cat = ifelse(df_full$cont_var == "LatinAmerica", 1,
                  ifelse(df_full$cont_var == "Africa", 2, 3))
# create match data.frame 
match_df = data.frame("uni_ID" = df_full$uni_ID,
                      "biome_car_cat" = biome_three_var,
                      "cont_cat" = cont_cat, 
                      stringsAsFactors = FALSE)
# merge the dataset
carb_stock_val = left_join(match_df, carbon_stock_dec, by = c("biome_car_cat" = "biome_category",
                                                              "cont_cat" = "contin_cat"))

carb_stock_val = carb_stock_val %>% select(uni_ID, biome_car_cat)

# merge the datasets
df_lat = inner_join(df_lat, carb_stock_val, by = "uni_ID")
df_afr = inner_join(df_afr, carb_stock_val, by = "uni_ID")
df_asi = inner_join(df_asi, carb_stock_val, by = "uni_ID")

df_lat$biome_num_f = factor(df_lat$biome_num)
df_afr$biome_num_f = factor(df_afr$biome_num)
df_asi$biome_num_f = factor(df_asi$biome_num)


# africa loss
reg_loss_afr <- glm(per_for_loss ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                           per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                           elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                         family="poisson", 
                         data=df_afr[(df_afr$per_forest_2000 > 0) & 
                                        (df_afr$biome_car_cat %in% c(1,2,3)), ])

# africa loss
reg_loss_asi <- glm(per_for_loss ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                      per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                      elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                    family="poisson", 
                    data=df_asi[(df_asi$per_forest_2000 > 0) & 
                                  (df_asi$biome_car_cat %in% c(1,2,3)), ])


# africa loss
reg_loss_lat <- glm(per_for_loss ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                      per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                      elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                    family="poisson", 
                    data=df_lat[(df_lat$per_forest_2000 > 0) & 
                                  (df_lat$biome_car_cat %in% c(1,2,3)), ])


rsq_vec = c(rsq(reg_loss_afr), rsq(reg_loss_asi), rsq(reg_loss_lat))
rsq_vec = round(rsq_vec, digits = 2)



### create stargazer results 
library(stargazer)
out_path = "D:/Research/global_refor_mac/results/regression_html/"

stargazer(reg_loss_afr, reg_loss_asi, reg_loss_lat,
          type = "html",
          column.labels = c("Africa", "Asia", "Latin America"), 
          title = "Continental land use change model of deforestation",
          out = paste0(out_path, "cont_land_use_cont_biome_excl_biomes.html"),
          p.auto = TRUE,
          align=FALSE, 
          digits = 5,
          font.size = "small",
          star.cutoffs=c(0.10,0.05,0.01),
          dep.var.labels=c("", "", ""),
          model.numbers = TRUE,
          multicolumn = FALSE,  
          perl = TRUE,
          keep = c(names(reg_loss_afr$coefficients)[c(2:11, 1)]),
          covariate.labels= c("Agricultural revenue ( hec)",
                              "Percent forest cover",
                              "Percent forest cover 2",
                              "Percent forest cover 3",
                              "Percent forest cover 4",
                              "Distance to city (km)",
                              "Slope (degree)",
                              "Elevation (m)",
                              "Multiple protected area (percent of cell)",
                              "Strict protected area (percent of cell)",
                              "Constant"),
          no.space = TRUE,
          keep.stat=c("n", "rsq"),
          single.row = FALSE,
          add.lines = list(c("R-squared", rsq_vec),
                           c("Biome FE", "Yes", "Yes", "Yes")))


reg_gain_afr = rxGlm(per_for_gain ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                        per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                        elev + iucnlow_2000 + iucnhigh_2000  + F(biome_num_f), 
                      data = df_afr[(df_afr$per_forest_2000 < 1) & 
                                      (df_afr$biome_car_cat %in% c(1,2,3)), ], 
                      family = poisson(link = log), 
                      coefLabelStyle = "R", 
                      maxIterations = 1000)

reg_gain_lat = rxGlm(per_for_gain ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                       per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                       elev + iucnlow_2000 + iucnhigh_2000  + F(biome_num_f), 
                     data = df_lat[(df_lat$per_forest_2000 < 1) & 
                                     (df_lat$biome_car_cat %in% c(1,2,3)), ], 
                     family = poisson(link = log), 
                     coefLabelStyle = "R", 
                     maxIterations = 1000)

reg_gain_asi = rxGlm(per_for_gain ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                       per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                       elev + iucnlow_2000 + iucnhigh_2000  + F(biome_num_f), 
                     data = df_asi[(df_asi$per_forest_2000 < 1) & 
                                     (df_asi$biome_car_cat %in% c(1,2,3)), ], 
                     family = poisson(link = log), 
                     coefLabelStyle = "R", 
                     maxIterations = 1000)

# save out the rds files 
saveRDS(reg_gain_afr, paste0("D:/Research/global_refor_mac/results/rds_reg_files/", "reg_gain_cont_v1_afr.rds"))
saveRDS(reg_gain_lat, paste0("D:/Research/global_refor_mac/results/rds_reg_files/", "reg_gain_cont_v1_lat.rds"))
saveRDS(reg_gain_asi, paste0("D:/Research/global_refor_mac/results/rds_reg_files/", "reg_gain_cont_v1_asi.rds"))



# africa loss
reg_gain_afr <- glm(per_for_gain ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                      per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                      elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                    family="poisson", 
                    data=df_afr[(df_afr$per_forest_2000 < 1) & 
                                  (df_afr$biome_car_cat %in% c(1,2,3)), ])

# africa loss
reg_gain_asi <- glm(per_for_gain ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                      per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                      elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                    family="poisson", 
                    data=df_asi[(df_asi$per_forest_2000 < 1) & 
                                  (df_asi$biome_car_cat %in% c(1,2,3)), ])


# africa loss
reg_gain_lat <- glm(per_for_gain ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                      per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                      elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                    family="poisson", 
                    data=df_lat[(df_lat$per_forest_2000 < 1) & 
                                  (df_lat$biome_car_cat %in% c(1,2,3)), ])

rsq_vec = c(rsq(reg_loss_afr), rsq(reg_loss_asi), rsq(reg_loss_lat))
rsq_vec = round(rsq_vec, digits = 2)


stargazer(reg_gain_afr, reg_gain_asi, reg_gain_lat,
          type = "html",
          column.labels = c("Africa", "Asia", "Latin America"), 
          title = "Continental land use change model of reforestation",
          out = paste0(out_path, "cont_reforestation_cont_biome_excl_biomes.html"),
          p.auto = TRUE,
          align=FALSE, 
          digits = 5,
          font.size = "small",
          star.cutoffs=c(0.10,0.05,0.01),
          dep.var.labels=c("", "", ""),
          model.numbers = TRUE,
          multicolumn = FALSE,  
          perl = TRUE,
          keep = c(names(reg_loss_afr$coefficients)[c(2:11, 1)]),
          covariate.labels= c("Agricultural revenue ( hec)",
                              "Percent forest cover",
                              "Percent forest cover 2",
                              "Percent forest cover 3",
                              "Percent forest cover 4",
                              "Distance to city (km)",
                              "Slope (degree)",
                              "Elevation (m)",
                              "Multiple protected area (percent of cell)",
                              "Strict protected area (percent of cell)",
                              "Constant"),
          no.space = TRUE,
          keep.stat=c("n", "rsq"),
          single.row = FALSE,
          add.lines = list(c("R-squared", rsq_vec),
                           c("Biome FE", "Yes", "Yes", "Yes")))


#####################################################################
### subtract oil plan reforestation 
df_plant2 = readRDS("D:/Research/global_refor_mac/tabulated_data/full_data/df_for_by_all_30_per.rds")
df_full = left_join(df_full, df_plant2, by = "uni_ID")
df_full2 = df_full

df_full2$per_for_gain = df_full2$per_for_gain - df_full2$per_gain_oil_palm
df_full2$per_for_loss = df_full2$per_for_loss - df_full2$per_loss_oil_palm

reg_loss_cont = rxGlm(per_for_loss ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                        per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                        elev + iucnlow_2000 + iucnhigh_2000 + F(biome_num) + F(cont_num_f), 
                      data = df_full2[(df_full2$per_forest_2000 > 0), ], 
                      family = poisson(link = log), 
                      coefLabelStyle = "R", 
                      maxIterations = 1000)

reg_gain_cont = rxGlm(per_for_gain ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                        per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                        elev + iucnlow_2000 + iucnhigh_2000  + F(biome_num) + F(cont_num_f), 
                      data = df_full2[df_full2$per_forest_2000 < 1, ], 
                      family = poisson(link = log), 
                      coefLabelStyle = "R", 
                      maxIterations = 1000)

saveRDS(reg_loss_cont, paste0("D:/Research/global_refor_mac/results/rds_reg_files/", "reg_loss_cont_oil_palm_excl.rds"))
saveRDS(reg_gain_cont, paste0("D:/Research/global_refor_mac/results/rds_reg_files/", "reg_gain_cont_oil_palm_excl.rds"))
rm(df_full2)
gc()

reg_gain_cont_glm <- glm(per_for_gain_clean ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                           per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                           elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f + cont_num_f, 
                         family="poisson", 
                         data=df_full[df_full$per_forest_2000 > 0, ])

library(stargazer)
out_path = "D:/Research/global_refor_mac/results/regression_html/"

stargazer(reg_gain_cont_glm,
          type = "html",
          title = "Land use change of reforestation (exclude oil palm)",
          out = paste0(out_path, "global_reforestation_cont_biome_excl_palmoil.html"),
          p.auto = TRUE,
          align=FALSE, 
          digits = 5,
          font.size = "small",
          star.cutoffs=c(0.10,0.05,0.01),
          dep.var.labels=c("Percent deforestation", "Percent reforestation"),
          model.numbers = TRUE,
          multicolumn = FALSE,  
          perl = TRUE,
          keep = c(names(reg_gain_cont_glm$coefficients)[c(2:11, 23,24,1)]),
          covariate.labels= c("Agricultural revenue ( hec)",
                              "Percent forest cover",
                              "Percent forest cover^2",
                              "Percent forest cover^3",
                              "Percent forest cover^4",
                              "Distance to city (km)",
                              "Slope (degree)",
                              "Elevation (m)",
                              "Multiple protected area (percent of cell)",
                              "Strict protected area (percent of cell)",
                              "Africa FE",
                              "Asia FE",
                              "Intercept"),
          no.space = TRUE,
          keep.stat=c("n", "rsq"),
          single.row = FALSE,
          add.lines = list(c("Biome FE", "Yes")))


# regresssion of palm oil 
reg_loss_cont_glm <- glm(per_for_loss_clean ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                           per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                           elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f + cont_num_f, 
                         family="poisson", 
                         data=df_full[df_full$per_forest_2000 > 0, ])


stargazer(reg_loss_cont_glm,
          type = "html",
          title = "Land use change of deforestation (exclude oil palm)",
          out = paste0(out_path, "global_deforestation_cont_biome_excl_palmoil.html"),
          p.auto = TRUE,
          align=FALSE, 
          digits = 5,
          font.size = "small",
          star.cutoffs=c(0.10,0.05,0.01),
          dep.var.labels=c("Percent deforestation", "Percent reforestation"),
          model.numbers = TRUE,
          multicolumn = FALSE,  
          perl = TRUE,
          keep = c(names(reg_gain_cont_glm$coefficients)[c(2:11, 23,24,1)]),
          covariate.labels= c("Agricultural revenue ( hec)",
                              "Percent forest cover",
                              "Percent forest cover^2",
                              "Percent forest cover^3",
                              "Percent forest cover^4",
                              "Distance to city (km)",
                              "Slope (degree)",
                              "Elevation (m)",
                              "Multiple protected area (percent of cell)",
                              "Strict protected area (percent of cell)",
                              "Africa FE",
                              "Asia FE",
                              "Intercept"),
          no.space = TRUE,
          keep.stat=c("n", "rsq"),
          single.row = FALSE,
          add.lines = list(c("Biome FE", "Yes")))


asi_reg_loss_cont_glm <- glm(per_for_loss_clean ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                           per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                           elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                         family="poisson", 
                         data=df_full[(df_full$per_forest_2000 > 0) & (df_full$cont_num_f == 3), ])


all_asi_reg_loss_cont_glm <- glm(per_for_loss ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                               per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                               elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                             family="poisson", 
                             data=df_full[(df_full$per_forest_2000 > 0) & (df_full$cont_num_f == 3), ])

stargazer(asi_reg_loss_cont_glm, all_asi_reg_loss_cont_glm,
          type = "html",
          title = "Land use change of deforestation (exclude oil palm)",
          out = paste0(out_path, "asi_deforestation_cont_biome_excl_palmoil.html"),
          p.auto = TRUE,
          align=FALSE, 
          column.labels = c("Exclude oil palm", "Standard"), 
          digits = 5,
          font.size = "small",
          star.cutoffs=c(0.10,0.05,0.01),
          dep.var.labels=c("Percent deforestation", "Percent reforestation"),
          model.numbers = TRUE,
          multicolumn = FALSE,  
          perl = TRUE,
          omit.stat = "all",
          keep = c(names(asi_reg_loss_cont_glm$coefficients)[c(2:11, 1)])),
          covariate.labels= c("Agricultural revenue ( hec)",
                              "Percent forest cover",
                              "Percent forest cover^2",
                              "Percent forest cover^3",
                              "Percent forest cover^4",
                              "Distance to city (km)",
                              "Slope (degree)",
                              "Elevation (m)",
                              "Multiple protected area (percent of cell)",
                              "Strict protected area (percent of cell)",
                              "Intercept")),
          no.space = TRUE,
          single.row = FALSE,
          add.lines = list(c("Biome FE", "Yes", "Yes")))



#########################################################################




#########################################################################

##### create country datasets
df_seven_cntr = df_full[df_full$country %in% c("Brazil",
                                               "Cambodia", 
                                               "Colombia", 
                                               "Indonesia", 
                                               "Liberia", 
                                               "Malaysia", 
                                               "Peru"), ]

df_seven_cntr$per_gain_timber = df_seven_cntr$per_gain_teak + df_seven_cntr$per_gain_softwood
df_seven_cntr$per_loss_timber = df_seven_cntr$per_loss_teak + df_seven_cntr$per_loss_softwood

reg_gain_oil_palm <- glm(per_gain_oil_palm ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                           per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                           elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                         family="poisson", 
                         data=df_seven_cntr[df_seven_cntr$per_forest_2000 < 1, ])

reg_gain_timber <- glm(per_gain_timber ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                           per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                           elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                         family="poisson", 
                         data=df_seven_cntr[df_seven_cntr$per_forest_2000 < 1, ])

reg_gain_other <- glm(per_gain_other ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                           per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                           elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                         family="poisson", 
                         data=df_seven_cntr[df_seven_cntr$per_forest_2000 < 1, ])

reg_gain_no_plant <- glm(per_gain_no_plant ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                           per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                           elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                         family="poisson", 
                         data=df_seven_cntr[df_seven_cntr$per_forest_2000 < 1, ])

reg_gain_all <- glm(per_for_gain ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                           per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                           elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                         family="poisson", 
                         data=df_seven_cntr[df_seven_cntr$per_forest_2000  < 1, ])


stargazer(reg_gain_oil_palm, reg_gain_timber, reg_gain_other, reg_gain_no_plant, reg_gain_all,
          type = "html",
          column.labels = c("Oil palm", "Timber", "Other", "No plantation", "All"), 
          title = "Seven countries land use change model of reforestation",
          out = paste0(out_path, "sev_cntry_reforestation_cont_biome_excl_biomes.html"),
          p.auto = TRUE,
          align=FALSE, 
          digits = 5,
          font.size = "small",
          star.cutoffs=c(0.10,0.05,0.01),
          dep.var.labels=c("", "", "", "", ""),
          model.numbers = TRUE,
          multicolumn = FALSE,  
          perl = TRUE,
          keep = c(names(reg_gain_oil_palm$coefficients)[c(2:11, 1)]),
          covariate.labels= c("Agricultural revenue ( hec)",
                              "Percent forest cover",
                              "Percent forest cover 2",
                              "Percent forest cover 3",
                              "Percent forest cover 4",
                              "Distance to city (km)",
                              "Slope (degree)",
                              "Elevation (m)",
                              "Multiple protected area (percent of cell)",
                              "Strict protected area (percent of cell)",
                              "Constant"),
          no.space = TRUE,
          keep.stat=c("n", "rsq"),
          single.row = FALSE,
          add.lines = list(c("Biome FE", "Yes", "Yes", "Yes", "Yes")))



reg_loss_oil_palm <- glm(per_loss_oil_palm ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                           per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                           elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                         family="poisson", 
                         data=df_seven_cntr[df_seven_cntr$per_forest_2000 > 0, ])

reg_loss_timber <- glm(per_loss_timber ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                         per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                         elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                       family="poisson", 
                       data=df_seven_cntr[df_seven_cntr$per_forest_2000 > 0, ])

reg_loss_other <- glm(per_loss_other ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                        per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                        elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                      family="poisson", 
                      data=df_seven_cntr[df_seven_cntr$per_forest_2000 > 0, ])

reg_loss_no_plant <- glm(per_loss_no_plant ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                           per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                           elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                         family="poisson", 
                         data=df_seven_cntr[df_seven_cntr$per_forest_2000 > 0, ])

reg_loss_all <- glm(per_for_loss ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                           per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                           elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                         family="poisson", 
                         data=df_seven_cntr[df_seven_cntr$per_forest_2000 > 0, ])


stargazer(reg_loss_oil_palm, reg_loss_timber, reg_loss_other, reg_loss_no_plant, reg_loss_all,
          type = "html",
          column.labels = c("Oil palm", "Timber", "Other", "No plantation", "All"), 
          title = "Seven countries land use change model of deforestation",
          out = paste0(out_path, "sev_cntry_deforestation_cont_biome_excl_biomes.html"),
          p.auto = TRUE,
          align=FALSE, 
          digits = 5,
          font.size = "small",
          star.cutoffs=c(0.10,0.05,0.01),
          dep.var.labels=c("", "", "", "", ""),
          model.numbers = TRUE,
          multicolumn = FALSE,  
          perl = TRUE,
          keep = c(names(reg_gain_oil_palm$coefficients)[c(2:11, 1)]),
          covariate.labels= c("Agricultural revenue ( hec)",
                              "Percent forest cover",
                              "Percent forest cover 2",
                              "Percent forest cover 3",
                              "Percent forest cover 4",
                              "Distance to city (km)",
                              "Slope (degree)",
                              "Elevation (m)",
                              "Multiple protected area (percent of cell)",
                              "Strict protected area (percent of cell)",
                              "Constant"),
          no.space = TRUE,
          keep.stat=c("n", "rsq"),
          single.row = FALSE,
          add.lines = list(c("Biome FE", "Yes", "Yes", "Yes", "Yes")))


##############################################################
### Brazil regression 

##### create country datasets
df_brazil = df_seven_cntr[df_seven_cntr$country %in% c("Brazil"), ]

reg_gain_oil_palm <- glm(per_gain_oil_palm ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                           per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                           elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                         family="poisson", 
                         data=df_brazil[df_brazil$per_forest_2000 < 1, ])

reg_gain_timber <- glm(per_gain_timber ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                         per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                         elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                       family="poisson", 
                       data=df_brazil[df_brazil$per_forest_2000 < 1, ])

reg_gain_other <- glm(per_gain_other ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                        per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                        elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                      family="poisson", 
                      data=df_brazil[df_brazil$per_forest_2000 < 1, ])

reg_gain_no_plant <- glm(per_gain_no_plant ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                           per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                           elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                         family="poisson", 
                         data=df_brazil[df_brazil$per_forest_2000 < 1, ])

reg_gain_all <- glm(per_for_gain ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                      per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                      elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                    family="poisson", 
                    data=df_brazil[df_brazil$per_forest_2000  < 1,  ])

# save out the rds files 
saveRDS(reg_gain_no_plant, paste0("D:/Research/global_refor_mac/results/rds_reg_files/", "reg_gain_nat_v1_brazil.rds"))
saveRDS(reg_gain_all, paste0("D:/Research/global_refor_mac/results/rds_reg_files/", "reg_gain_all_v1_brazil.rds"))

rsq_vec = c(rsq(reg_gain_oil_palm), rsq(reg_gain_timber), rsq(reg_gain_other),
            rsq(reg_gain_no_plant), rsq(reg_gain_all))
rsq_vec = round(rsq_vec, digits = 2)



stargazer(reg_gain_oil_palm, reg_gain_timber, reg_gain_other, reg_gain_no_plant, reg_gain_all,
          type = "html",
          column.labels = c("Oil palm", "Timber", "Other", "No plantation", "All"), 
          title = "Brazil land use change model of reforestation",
          out = paste0(out_path, "brazil_reforestation_cont_biome_excl_biomes.html"),
          p.auto = TRUE,
          align=FALSE, 
          digits = 5,
          font.size = "small",
          star.cutoffs=c(0.10,0.05,0.01),
          dep.var.labels=c("", "", "", "", ""),
          model.numbers = TRUE,
          multicolumn = FALSE,  
          perl = TRUE,
          keep = c(names(reg_gain_oil_palm$coefficients)[c(2:11, 1)]),
          covariate.labels= c("Agricultural revenue ( hec)",
                              "Percent forest cover",
                              "Percent forest cover 2",
                              "Percent forest cover 3",
                              "Percent forest cover 4",
                              "Distance to city (km)",
                              "Slope (degree)",
                              "Elevation (m)",
                              "Multiple protected area (percent of cell)",
                              "Strict protected area (percent of cell)",
                              "Constant"),
          no.space = TRUE,
          keep.stat=c("n", "rsq"),
          single.row = FALSE,
          add.lines = list(c("R-squared", rsq_vec),
                           c("Biome FE", "Yes", "Yes", "Yes", "Yes", "Yes")))



reg_loss_oil_palm <- glm(per_loss_oil_palm ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                           per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                           elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                         family="poisson", 
                         data=df_brazil[df_brazil$per_forest_2000 > 0, ])

reg_loss_timber <- glm(per_loss_timber ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                         per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                         elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                       family="poisson", 
                       data=df_brazil[df_brazil$per_forest_2000 > 0, ])

reg_loss_other <- glm(per_loss_other ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                        per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                        elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                      family="poisson", 
                      data=df_brazil[df_brazil$per_forest_2000 > 0, ])

reg_loss_no_plant <- glm(per_loss_no_plant ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                           per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                           elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                         family="poisson", 
                         data=df_brazil[df_brazil$per_forest_2000 > 0, ])

reg_loss_all <- glm(per_for_loss ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                      per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                      elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                    family="poisson", 
                    data=df_brazil[df_brazil$per_forest_2000 > 0, ])


stargazer(reg_loss_oil_palm, reg_loss_timber, reg_loss_other, reg_loss_no_plant, reg_loss_all,
          type = "html",
          column.labels = c("Oil palm", "Timber", "Other", "No plantation", "All"), 
          title = "Brazil land use change model of deforestation",
          out = paste0(out_path, "brazil_deforestation_cont_biome_excl_biomes.html"),
          p.auto = TRUE,
          align=FALSE, 
          digits = 5,
          font.size = "small",
          star.cutoffs=c(0.10,0.05,0.01),
          dep.var.labels=c("", "", "", "", ""),
          model.numbers = TRUE,
          multicolumn = FALSE,  
          perl = TRUE,
          keep = c(names(reg_gain_oil_palm$coefficients)[c(2:11, 1)]),
          covariate.labels= c("Agricultural revenue ( hec)",
                              "Percent forest cover",
                              "Percent forest cover 2",
                              "Percent forest cover 3",
                              "Percent forest cover 4",
                              "Distance to city (km)",
                              "Slope (degree)",
                              "Elevation (m)",
                              "Multiple protected area (percent of cell)",
                              "Strict protected area (percent of cell)",
                              "Constant"),
          no.space = TRUE,
          keep.stat=c("n", "rsq"),
          single.row = FALSE,
          add.lines = list(c("Biome FE", "Yes", "Yes", "Yes", "Yes")))



##############################################################
### Cambodia regression 

##### create country datasets
df_camb = df_seven_cntr[df_seven_cntr$country %in% c("Cambodia"), ]

reg_gain_oil_palm <- glm(per_gain_oil_palm ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                           per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                           elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                         family="poisson", 
                         data=df_camb[df_camb$per_forest_2000  < 1, ])

reg_gain_timber <- glm(per_gain_timber ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                         per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                         elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                       family="poisson", 
                       data=df_camb[df_camb$per_forest_2000  < 1, ])

reg_gain_other <- glm(per_gain_other ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                        per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                        elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                      family="poisson", 
                      data=df_camb[df_camb$per_forest_2000  < 1, ])

reg_gain_no_plant <- glm(per_gain_no_plant ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                           per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                           elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                         family="poisson", 
                         data=df_camb[df_camb$per_forest_2000  < 1, ])

reg_gain_all <- glm(per_for_gain ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                      per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                      elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                    family="poisson", 
                    data=df_camb[df_camb$per_forest_2000 < 1, ])


stargazer(reg_gain_oil_palm, reg_gain_timber, reg_gain_other, reg_gain_no_plant, reg_gain_all,
          type = "html",
          column.labels = c("Oil palm", "Timber", "Other", "No plantation", "All"), 
          title = "Cambodia land use change model of reforestation",
          out = paste0(out_path, "cambodia_reforestation_cont_biome_excl_biomes.html"),
          p.auto = TRUE,
          align=FALSE, 
          digits = 5,
          font.size = "small",
          star.cutoffs=c(0.10,0.05,0.01),
          dep.var.labels=c("", "", "", "", ""),
          model.numbers = TRUE,
          multicolumn = FALSE,  
          perl = TRUE,
          keep = c(names(reg_gain_oil_palm$coefficients)[c(2:11, 1)]),
          covariate.labels= c("Agricultural revenue ( hec)",
                              "Percent forest cover",
                              "Percent forest cover 2",
                              "Percent forest cover 3",
                              "Percent forest cover 4",
                              "Distance to city (km)",
                              "Slope (degree)",
                              "Elevation (m)",
                              "Multiple protected area (percent of cell)",
                              "Strict protected area (percent of cell)",
                              "Constant"),
          no.space = TRUE,
          keep.stat=c("n", "rsq"),
          single.row = FALSE,
          add.lines = list(c("Biome FE", "Yes", "Yes", "Yes", "Yes")))



reg_loss_oil_palm <- glm(per_loss_oil_palm ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                           per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                           elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                         family="poisson", 
                         data=df_camb[df_camb$per_forest_2000 > 0, ])

reg_loss_timber <- glm(per_loss_timber ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                         per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                         elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                       family="poisson", 
                       data=df_camb[df_camb$per_forest_2000 > 0, ])

reg_loss_other <- glm(per_loss_other ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                        per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                        elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                      family="poisson", 
                      data=df_camb[df_camb$per_forest_2000 > 0, ])

reg_loss_no_plant <- glm(per_loss_no_plant ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                           per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                           elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                         family="poisson", 
                         data=df_camb[df_camb$per_forest_2000 > 0, ])

reg_loss_all <- glm(per_for_loss ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                      per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                      elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                    family="poisson", 
                    data=df_camb[df_camb$per_forest_2000 > 0, ])


stargazer(reg_loss_oil_palm, reg_loss_timber, reg_loss_other, reg_loss_no_plant, reg_loss_all,
          type = "html",
          column.labels = c("Oil palm", "Timber", "Other", "No plantation", "All"), 
          title = "Cambodia land use change model of deforestation",
          out = paste0(out_path, "cambodia_deforestation_cont_biome_excl_biomes.html"),
          p.auto = TRUE,
          align=FALSE, 
          digits = 5,
          font.size = "small",
          star.cutoffs=c(0.10,0.05,0.01),
          dep.var.labels=c("", "", "", "", ""),
          model.numbers = TRUE,
          multicolumn = FALSE,  
          perl = TRUE,
          keep = c(names(reg_gain_oil_palm$coefficients)[c(2:11, 1)]),
          covariate.labels= c("Agricultural revenue ( hec)",
                              "Percent forest cover",
                              "Percent forest cover 2",
                              "Percent forest cover 3",
                              "Percent forest cover 4",
                              "Distance to city (km)",
                              "Slope (degree)",
                              "Elevation (m)",
                              "Multiple protected area (percent of cell)",
                              "Strict protected area (percent of cell)",
                              "Constant"),
          no.space = TRUE,
          keep.stat=c("n", "rsq"),
          single.row = FALSE,
          add.lines = list(c("Biome FE", "Yes", "Yes", "Yes", "Yes")))












##############################################################
### Cambodia regression 

##### create country datasets
df_colomb = df_seven_cntr[df_seven_cntr$country %in% c("Colombia"), ]

reg_gain_oil_palm <- glm(per_gain_oil_palm ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                           per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                           elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                         family="poisson", 
                         data=df_colomb[df_colomb$per_forest_2000  < 1, ])

reg_gain_timber <- glm(per_gain_timber ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                         per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                         elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                       family="poisson", 
                       data=df_colomb[df_colomb$per_forest_2000 < 1, ])

reg_gain_other <- glm(per_gain_other ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                        per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                        elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                      family="poisson", 
                      data=df_colomb[df_colomb$per_forest_2000  < 1, ])

reg_gain_no_plant <- glm(per_gain_no_plant ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                           per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                           elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                         family="poisson", 
                         data=df_colomb[df_colomb$per_forest_2000  < 1, ])

reg_gain_all <- glm(per_for_gain ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                      per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                      elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                    family="poisson", 
                    data=df_colomb[df_colomb$per_forest_2000  < 1, ])


stargazer(reg_gain_oil_palm, reg_gain_timber, reg_gain_other, reg_gain_no_plant, reg_gain_all,
          type = "html",
          column.labels = c("Oil palm", "Timber", "Other", "No plantation", "All"), 
          title = "Colombia land use change model of reforestation",
          out = paste0(out_path, "colombia_reforestation_cont_biome_excl_biomes.html"),
          p.auto = TRUE,
          align=FALSE, 
          digits = 5,
          font.size = "small",
          star.cutoffs=c(0.10,0.05,0.01),
          dep.var.labels=c("", "", "", "", ""),
          model.numbers = TRUE,
          multicolumn = FALSE,  
          perl = TRUE,
          keep = c(names(reg_gain_oil_palm$coefficients)[c(2:11, 1)]),
          covariate.labels= c("Agricultural revenue ( hec)",
                              "Percent forest cover",
                              "Percent forest cover 2",
                              "Percent forest cover 3",
                              "Percent forest cover 4",
                              "Distance to city (km)",
                              "Slope (degree)",
                              "Elevation (m)",
                              "Multiple protected area (percent of cell)",
                              "Strict protected area (percent of cell)",
                              "Constant"),
          no.space = TRUE,
          keep.stat=c("n", "rsq"),
          single.row = FALSE,
          add.lines = list(c("Biome FE", "Yes", "Yes", "Yes", "Yes")))



reg_loss_oil_palm <- glm(per_loss_oil_palm ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                           per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                           elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                         family="poisson", 
                         data=df_colomb[df_colomb$per_forest_2000 > 0, ])

reg_loss_timber <- glm(per_loss_timber ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                         per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                         elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                       family="poisson", 
                       data=df_colomb[df_colomb$per_forest_2000 > 0, ])

reg_loss_other <- glm(per_loss_other ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                        per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                        elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                      family="poisson", 
                      data=df_colomb[df_colomb$per_forest_2000 > 0, ])

reg_loss_no_plant <- glm(per_loss_no_plant ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                           per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                           elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                         family="poisson", 
                         data=df_colomb[df_colomb$per_forest_2000 > 0, ])

reg_loss_all <- glm(per_for_loss ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                      per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                      elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                    family="poisson", 
                    data=df_colomb[df_colomb$per_forest_2000 > 0, ])


stargazer(reg_loss_oil_palm, reg_loss_timber, reg_loss_other, reg_loss_no_plant, reg_loss_all,
          type = "html",
          column.labels = c("Oil palm", "Timber", "Other", "No plantation", "All"), 
          title = "Colombia land use change model of deforestation",
          out = paste0(out_path, "colombia_deforestation_cont_biome_excl_biomes.html"),
          p.auto = TRUE,
          align=FALSE, 
          digits = 5,
          font.size = "small",
          star.cutoffs=c(0.10,0.05,0.01),
          dep.var.labels=c("", "", "", "", ""),
          model.numbers = TRUE,
          multicolumn = FALSE,  
          perl = TRUE,
          keep = c(names(reg_gain_oil_palm$coefficients)[c(2:11, 1)]),
          covariate.labels= c("Agricultural revenue ( hec)",
                              "Percent forest cover",
                              "Percent forest cover 2",
                              "Percent forest cover 3",
                              "Percent forest cover 4",
                              "Distance to city (km)",
                              "Slope (degree)",
                              "Elevation (m)",
                              "Multiple protected area (percent of cell)",
                              "Strict protected area (percent of cell)",
                              "Constant"),
          no.space = TRUE,
          keep.stat=c("n", "rsq"),
          single.row = FALSE,
          add.lines = list(c("Biome FE", "Yes", "Yes", "Yes", "Yes")))




##############################################################
### Cambodia regression 

##### create country datasets
df_idn = df_seven_cntr[df_seven_cntr$country %in% c("Indonesia"), ]

reg_gain_oil_palm <- glm(per_gain_oil_palm ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                           per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                           elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                         family="poisson", 
                         data=df_idn[df_idn$per_forest_2000  < 1, ])

reg_gain_timber <- glm(per_gain_timber ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                         per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                         elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                       family="poisson", 
                       data=df_idn[df_idn$per_forest_2000  < 1, ])

reg_gain_other <- glm(per_gain_other ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                        per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                        elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                      family="poisson", 
                      data=df_idn[df_idn$per_forest_2000  < 1, ])

reg_gain_no_plant <- glm(per_gain_no_plant ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                           per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                           elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                         family="poisson", 
                         data=df_idn[df_idn$per_forest_2000  < 1, ])

reg_gain_all <- glm(per_for_gain ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                      per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                      elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                    family="poisson", 
                    data=df_idn[df_idn$per_forest_2000  < 1, ])

rsq_vec = c(rsq(reg_gain_oil_palm), rsq(reg_gain_timber), rsq(reg_gain_other),
            rsq(reg_gain_no_plant), rsq(reg_gain_all))
rsq_vec = round(rsq_vec, digits = 2)


stargazer(reg_gain_oil_palm, reg_gain_timber, reg_gain_other, reg_gain_no_plant, reg_gain_all,
          type = "html",
          column.labels = c("Oil palm", "Timber", "Other", "No plantation", "All"), 
          title = "Indonesia land use change model of reforestation",
          out = paste0(out_path, "indonesia_reforestation_cont_biome_excl_biomes.html"),
          p.auto = TRUE,
          align=FALSE, 
          digits = 5,
          font.size = "small",
          star.cutoffs=c(0.10,0.05,0.01),
          dep.var.labels=c("", "", "", "", ""),
          model.numbers = TRUE,
          multicolumn = FALSE,  
          perl = TRUE,
          keep = c(names(reg_gain_oil_palm$coefficients)[c(2:11, 1)]),
          covariate.labels= c("Agricultural revenue ( hec)",
                              "Percent forest cover",
                              "Percent forest cover 2",
                              "Percent forest cover 3",
                              "Percent forest cover 4",
                              "Distance to city (km)",
                              "Slope (degree)",
                              "Elevation (m)",
                              "Multiple protected area (percent of cell)",
                              "Strict protected area (percent of cell)",
                              "Constant"),
          no.space = TRUE,
          keep.stat=c("n", "rsq"),
          single.row = FALSE,
          add.lines = list(c("R-squared", rsq_vec),
                           c("Biome FE", "Yes", "Yes", "Yes", "Yes")))



reg_loss_oil_palm <- glm(per_loss_oil_palm ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                           per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                           elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                         family="poisson", 
                         data=df_idn[df_idn$per_forest_2000 > 0, ])

reg_loss_timber <- glm(per_loss_timber ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                         per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                         elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                       family="poisson", 
                       data=df_idn[df_idn$per_forest_2000 > 0, ])

reg_loss_other <- glm(per_loss_other ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                        per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                        elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                      family="poisson", 
                      data=df_idn[df_idn$per_forest_2000 > 0, ])

reg_loss_no_plant <- glm(per_loss_no_plant ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                           per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                           elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                         family="poisson", 
                         data=df_idn[df_idn$per_forest_2000 > 0, ])

reg_loss_all <- glm(per_for_loss ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                      per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                      elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                    family="poisson", 
                    data=df_idn[df_idn$per_forest_2000 > 0, ])


stargazer(reg_loss_oil_palm, reg_loss_timber, reg_loss_other, reg_loss_no_plant, reg_loss_all,
          type = "html",
          column.labels = c("Oil palm", "Timber", "Other", "No plantation", "All"), 
          title = "Indonesia land use change model of deforestation",
          out = paste0(out_path, "indonesia_deforestation_cont_biome_excl_biomes.html"),
          p.auto = TRUE,
          align=FALSE, 
          digits = 5,
          font.size = "small",
          star.cutoffs=c(0.10,0.05,0.01),
          dep.var.labels=c("", "", "", "", ""),
          model.numbers = TRUE,
          multicolumn = FALSE,  
          perl = TRUE,
          keep = c(names(reg_gain_oil_palm$coefficients)[c(2:11, 1)]),
          covariate.labels= c("Agricultural revenue ( hec)",
                              "Percent forest cover",
                              "Percent forest cover 2",
                              "Percent forest cover 3",
                              "Percent forest cover 4",
                              "Distance to city (km)",
                              "Slope (degree)",
                              "Elevation (m)",
                              "Multiple protected area (percent of cell)",
                              "Strict protected area (percent of cell)",
                              "Constant"),
          no.space = TRUE,
          keep.stat=c("n", "rsq"),
          single.row = FALSE,
          add.lines = list(c("Biome FE", "Yes", "Yes", "Yes", "Yes")))





##############################################################
### Liberia regression 

##### create country datasets
df_lib = df_seven_cntr[df_seven_cntr$country %in% c("Liberia"), ]

reg_gain_oil_palm <- glm(per_gain_oil_palm ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                           per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                           elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                         family="poisson", 
                         data=df_lib[df_lib$per_forest_2000 < 1, ])

reg_gain_timber <- glm(per_gain_timber ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                         per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                         elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                       family="poisson", 
                       data=df_lib[df_lib$per_forest_2000  < 1, ])

reg_gain_other <- glm(per_gain_other ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                        per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                        elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                      family="poisson", 
                      data=df_lib[df_lib$per_forest_2000  < 1, ])

reg_gain_no_plant <- glm(per_gain_no_plant ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                           per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                           elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                         family="poisson", 
                         data=df_lib[df_lib$per_forest_2000  < 1, ])

reg_gain_all <- glm(per_for_gain ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                      per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                      elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                    family="poisson", 
                    data=df_lib[df_lib$per_forest_2000  < 1, ])


stargazer(reg_gain_oil_palm, reg_gain_other, reg_gain_no_plant, reg_gain_all,
          type = "html",
          column.labels = c("Oil palm", "Other", "No plantation", "All"), 
          title = "Liberia land use change model of reforestation",
          out = paste0(out_path, "liberia_reforestation_cont_biome_excl_biomes.html"),
          p.auto = TRUE,
          align=FALSE, 
          digits = 5,
          font.size = "small",
          star.cutoffs=c(0.10,0.05,0.01),
          dep.var.labels=c("", "", "", "", ""),
          model.numbers = TRUE,
          multicolumn = FALSE,  
          perl = TRUE,
          keep = c(names(reg_gain_oil_palm$coefficients)[c(2:11, 1)]),
          covariate.labels= c("Agricultural revenue ( hec)",
                              "Percent forest cover",
                              "Percent forest cover 2",
                              "Percent forest cover 3",
                              "Percent forest cover 4",
                              "Distance to city (km)",
                              "Slope (degree)",
                              "Elevation (m)",
                              "Multiple protected area (percent of cell)",
                              "Strict protected area (percent of cell)",
                              "Constant"),
          no.space = TRUE,
          keep.stat=c("n", "rsq"),
          single.row = FALSE,
          add.lines = list(c("Biome FE", "Yes", "Yes", "Yes", "Yes")))



reg_loss_oil_palm <- glm(per_loss_oil_palm ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                           per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                           elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                         family="poisson", 
                         data=df_lib[df_lib$per_forest_2000 > 0, ])

reg_loss_other <- glm(per_loss_other ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                        per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                        elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                      family="poisson", 
                      data=df_lib[df_lib$per_forest_2000 > 0, ])

reg_loss_no_plant <- glm(per_loss_no_plant ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                           per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                           elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                         family="poisson", 
                         data=df_lib[df_lib$per_forest_2000 > 0, ])

reg_loss_all <- glm(per_for_loss ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                      per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                      elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                    family="poisson", 
                    data=df_lib[df_lib$per_forest_2000 > 0, ])


stargazer(reg_loss_oil_palm, reg_loss_other, reg_loss_no_plant, reg_loss_all,
          type = "html",
          column.labels = c("Oil palm", "Timber", "Other", "No plantation", "All"), 
          title = "Liberia land use change model of deforestation",
          out = paste0(out_path, "liberia_deforestation_cont_biome_excl_biomes.html"),
          p.auto = TRUE,
          align=FALSE, 
          digits = 5,
          font.size = "small",
          star.cutoffs=c(0.10,0.05,0.01),
          dep.var.labels=c("", "", "", "", ""),
          model.numbers = TRUE,
          multicolumn = FALSE,  
          perl = TRUE,
          keep = c(names(reg_gain_oil_palm$coefficients)[c(2:11, 1)]),
          covariate.labels= c("Agricultural revenue ( hec)",
                              "Percent forest cover",
                              "Percent forest cover 2",
                              "Percent forest cover 3",
                              "Percent forest cover 4",
                              "Distance to city (km)",
                              "Slope (degree)",
                              "Elevation (m)",
                              "Multiple protected area (percent of cell)",
                              "Strict protected area (percent of cell)",
                              "Constant"),
          no.space = TRUE,
          keep.stat=c("n", "rsq"),
          single.row = FALSE,
          add.lines = list(c("Biome FE", "Yes", "Yes", "Yes", "Yes")))




##############################################################
### Cambodia regression 

##### create country datasets
df_mal = df_seven_cntr[df_seven_cntr$country %in% c("Malaysia"), ]

reg_gain_oil_palm <- glm(per_gain_oil_palm ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                           per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                           elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                         family="poisson", 
                         data=df_mal[df_mal$per_forest_2000 < 1, ])

reg_gain_timber <- glm(per_gain_timber ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                         per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                         elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                       family="poisson", 
                       data=df_mal[df_mal$per_forest_2000  < 1, ])

reg_gain_other <- glm(per_gain_other ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                        per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                        elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                      family="poisson", 
                      data=df_mal[df_mal$per_forest_2000  < 1, ])

reg_gain_no_plant <- glm(per_gain_no_plant ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                           per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                           elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                         family="poisson", 
                         data=df_mal[df_mal$per_forest_2000  < 1, ])

reg_gain_all <- glm(per_for_gain ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                      per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                      elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                    family="poisson", 
                    data=df_mal[df_mal$per_forest_2000  < 1, ])


stargazer(reg_gain_oil_palm, reg_gain_timber, reg_gain_other, reg_gain_no_plant, reg_gain_all,
          type = "html",
          column.labels = c("Oil palm", "Timber", "Other", "No plantation", "All"), 
          title = "Malaysia land use change model of reforestation",
          out = paste0(out_path, "malaysia_reforestation_cont_biome_excl_biomes.html"),
          p.auto = TRUE,
          align=FALSE, 
          digits = 5,
          font.size = "small",
          star.cutoffs=c(0.10,0.05,0.01),
          dep.var.labels=c("", "", "", "", ""),
          model.numbers = TRUE,
          multicolumn = FALSE,  
          perl = TRUE,
          keep = c(names(reg_gain_oil_palm$coefficients)[c(2:11, 1)]),
          covariate.labels= c("Agricultural revenue ( hec)",
                              "Percent forest cover",
                              "Percent forest cover 2",
                              "Percent forest cover 3",
                              "Percent forest cover 4",
                              "Distance to city (km)",
                              "Slope (degree)",
                              "Elevation (m)",
                              "Multiple protected area (percent of cell)",
                              "Strict protected area (percent of cell)",
                              "Constant"),
          no.space = TRUE,
          keep.stat=c("n", "rsq"),
          single.row = FALSE,
          add.lines = list(c("Biome FE", "Yes", "Yes", "Yes", "Yes")))



reg_loss_oil_palm <- glm(per_loss_oil_palm ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                           per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                           elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                         family="poisson", 
                         data=df_mal[df_mal$per_forest_2000 > 0, ])

reg_loss_timber <- glm(per_loss_timber ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                         per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                         elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                       family="poisson", 
                       data=df_mal[df_mal$per_forest_2000 > 0, ])

reg_loss_other <- glm(per_loss_other ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                        per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                        elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                      family="poisson", 
                      data=df_mal[df_mal$per_forest_2000 > 0, ])

reg_loss_no_plant <- glm(per_loss_no_plant ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                           per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                           elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                         family="poisson", 
                         data=df_mal[df_mal$per_forest_2000 > 0, ])

reg_loss_all <- glm(per_for_loss ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                      per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                      elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                    family="poisson", 
                    data=df_mal[df_mal$per_forest_2000 > 0, ])


stargazer(reg_loss_oil_palm, reg_loss_timber, reg_loss_other, reg_loss_no_plant, reg_loss_all,
          type = "html",
          column.labels = c("Oil palm", "Timber", "Other", "No plantation", "All"), 
          title = "Malaysia land use change model of deforestation",
          out = paste0(out_path, "malaysia_deforestation_cont_biome_excl_biomes.html"),
          p.auto = TRUE,
          align=FALSE, 
          digits = 5,
          font.size = "small",
          star.cutoffs=c(0.10,0.05,0.01),
          dep.var.labels=c("", "", "", "", ""),
          model.numbers = TRUE,
          multicolumn = FALSE,  
          perl = TRUE,
          keep = c(names(reg_gain_oil_palm$coefficients)[c(2:11, 1)]),
          covariate.labels= c("Agricultural revenue ( hec)",
                              "Percent forest cover",
                              "Percent forest cover 2",
                              "Percent forest cover 3",
                              "Percent forest cover 4",
                              "Distance to city (km)",
                              "Slope (degree)",
                              "Elevation (m)",
                              "Multiple protected area (percent of cell)",
                              "Strict protected area (percent of cell)",
                              "Constant"),
          no.space = TRUE,
          keep.stat=c("n", "rsq"),
          single.row = FALSE,
          add.lines = list(c("Biome FE", "Yes", "Yes", "Yes", "Yes")))



##############################################################
### Cambodia regression 

##### create country datasets
df_peru = df_seven_cntr[df_seven_cntr$country %in% c("Peru"), ]

reg_gain_oil_palm <- glm(per_gain_oil_palm ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                           per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                           elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                         family="poisson", 
                         data=df_peru[df_peru$per_forest_2000  < 1, ])

reg_gain_timber <- glm(per_gain_timber ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                         per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                         elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                       family="poisson", 
                       data=df_peru[df_peru$per_forest_2000  < 1, ])

reg_gain_other <- glm(per_gain_other ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                        per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                        elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                      family="poisson", 
                      data=df_peru[df_peru$per_forest_2000  < 1, ])

reg_gain_no_plant <- glm(per_gain_no_plant ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                           per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                           elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                         family="poisson", 
                         data=df_peru[df_peru$per_forest_2000  < 1, ])

reg_gain_all <- glm(per_for_gain ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                      per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                      elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                    family="poisson", 
                    data=df_peru[df_peru$per_forest_2000 < 1, ])


stargazer(reg_gain_oil_palm, reg_gain_timber, reg_gain_other, reg_gain_no_plant, reg_gain_all,
          type = "html",
          column.labels = c("Oil palm", "Timber", "Other", "No plantation", "All"), 
          title = "Peru land use change model of reforestation",
          out = paste0(out_path, "peru_reforestation_cont_biome_excl_biomes.html"),
          p.auto = TRUE,
          align=FALSE, 
          digits = 5,
          font.size = "small",
          star.cutoffs=c(0.10,0.05,0.01),
          dep.var.labels=c("", "", "", "", ""),
          model.numbers = TRUE,
          multicolumn = FALSE,  
          perl = TRUE,
          keep = c(names(reg_gain_oil_palm$coefficients)[c(2:11, 1)]),
          covariate.labels= c("Agricultural revenue ( hec)",
                              "Percent forest cover",
                              "Percent forest cover 2",
                              "Percent forest cover 3",
                              "Percent forest cover 4",
                              "Distance to city (km)",
                              "Slope (degree)",
                              "Elevation (m)",
                              "Multiple protected area (percent of cell)",
                              "Strict protected area (percent of cell)",
                              "Constant"),
          no.space = TRUE,
          keep.stat=c("n", "rsq"),
          single.row = FALSE,
          add.lines = list(c("Biome FE", "Yes", "Yes", "Yes", "Yes")))



reg_loss_oil_palm <- glm(per_loss_oil_palm ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                           per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                           elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                         family="poisson", 
                         data=df_peru[df_peru$per_forest_2000 > 0, ])

reg_loss_timber <- glm(per_loss_timber ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                         per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                         elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                       family="poisson", 
                       data=df_peru[df_peru$per_forest_2000 > 0, ])

reg_loss_other <- glm(per_loss_other ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                        per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                        elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                      family="poisson", 
                      data=df_peru[df_peru$per_forest_2000 > 0, ])

reg_loss_no_plant <- glm(per_loss_no_plant ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                           per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                           elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                         family="poisson", 
                         data=df_peru[df_peru$per_forest_2000 > 0, ])

reg_loss_all <- glm(per_for_loss ~ agrev_dec_max_00_10 + per_forest_2000 + per_forest_2000_2 +
                      per_forest_2000_3 + per_forest_2000_4 + dis_city_750k_near + slope +
                      elev + iucnlow_2000 + iucnhigh_2000 + biome_num_f, 
                    family="poisson", 
                    data=df_peru[df_peru$per_forest_2000 > 0, ])


stargazer(reg_loss_oil_palm, reg_loss_timber, reg_loss_other, reg_loss_no_plant, reg_loss_all,
          type = "html",
          column.labels = c("Oil palm", "Timber", "Other", "No plantation", "All"), 
          title = "Peru land use change model of deforestation",
          out = paste0(out_path, "peru_deforestation_cont_biome_excl_biomes.html"),
          p.auto = TRUE,
          align=FALSE, 
          digits = 5,
          font.size = "small",
          star.cutoffs=c(0.10,0.05,0.01),
          dep.var.labels=c("", "", "", "", ""),
          model.numbers = TRUE,
          multicolumn = FALSE,  
          perl = TRUE,
          keep = c(names(reg_gain_oil_palm$coefficients)[c(2:11, 1)]),
          covariate.labels= c("Agricultural revenue ( hec)",
                              "Percent forest cover",
                              "Percent forest cover 2",
                              "Percent forest cover 3",
                              "Percent forest cover 4",
                              "Distance to city (km)",
                              "Slope (degree)",
                              "Elevation (m)",
                              "Multiple protected area (percent of cell)",
                              "Strict protected area (percent of cell)",
                              "Constant"),
          no.space = TRUE,
          keep.stat=c("n", "rsq"),
          single.row = FALSE,
          add.lines = list(c("Biome FE", "Yes", "Yes", "Yes", "Yes")))

# create statistics by country 
df_out = df_seven_cntr %>% group_by(country) %>% mutate(gain_oil_palm = per_gain_oil_palm * shp_size,
                                               gain_timber = per_gain_timber * shp_size,
                                               gain_other = per_gain_other * shp_size,
                                               gain_no_plant = per_gain_no_plant * shp_size) %>%
                                        summarize(gain_oil_palm = sum(gain_oil_palm, na.rm = TRUE),
                                                  gain_timber = sum(gain_timber, na.rm = TRUE),
                                                  gain_other = sum(gain_other, na.rm = TRUE),
                                                  gain_no_plant = sum(gain_no_plant, na.rm = TRUE))


write.csv(df_out, paste0(out_path, "refor_by_cat_sev_cntry.csv"))


##################################################################################################
#### create reforestation percent gain and percent forest cover

# create statistics by country 
df_out = df_seven_cntr %>% group_by(country) %>% mutate(for_oil_palm = (per_forest_oil_palm + per_gain_oil_palm) * shp_size,
                                                        gain_oil_palm = per_gain_oil_palm * shp_size,
                                                        for_timber = (per_forest_teak + per_forest_softwood) * shp_size,
                                                        gain_timber = per_gain_timber * shp_size,
                                                        for_other = (per_forest_other + per_gain_other) * shp_size,
                                                        gain_other = per_gain_other * shp_size,
                                                        for_no_plant = (per_forest_no_plant + per_gain_no_plant) * shp_size,
                                                        gain_no_plant = per_gain_no_plant * shp_size) %>%
  summarize(for_oil_palm = sum(for_oil_palm, na.rm = TRUE),
            gain_oil_palm = sum(gain_oil_palm, na.rm = TRUE),
            for_timber = sum(for_timber, na.rm = TRUE),
            gain_timber = sum(gain_timber, na.rm = TRUE),
            for_other = sum(for_other, na.rm = TRUE),
            gain_other = sum(gain_other, na.rm = TRUE),
            for_no_plant = sum(for_no_plant, na.rm = TRUE),
            gain_no_plant = sum(gain_no_plant, na.rm = TRUE))


write.csv(df_out, paste0(out_path, "forest_cover_refor_by_cat_sev_cntry.csv"))
